Przyglądnijmy się teraz innej klasie problemów obliczeniowych, tak zwanemu problemowi konwekcji-dyfuzji. W problemie tym równanie różniczkowe w tak zwanej formie silnej dane jest wzorem
\( \beta_x\frac{\partial u}{\partial x }+\beta_y\frac{\partial u }{\partial y }-\epsilon \left(\frac{\partial^2 u }{\partial x^2 }+\frac{\partial^2 u }{\partial y^2 }\right)=f(x,y) \)
Szukamy pola koncentracji, np. zanieczyszczeń \( \Omega\in (x,y) \rightarrow u(x,y) \in R \) na obszarze \( \Omega \). Wartość \( u(x,y) \) oznacza koncentracje zanieczyszczeń w punkcie \( (x,y) \).
Człon \( \beta_x(x,y)\frac{\partial u }{\partial x }+\beta_y(x,y)\frac{\partial u }{\partial y } \) oznacza tutaj adwekcje, czyli roznoszenie zanieczyszczeń przez wiatr.
Składowe \( (\beta_x(x,y),\beta_y(x,y)) \) oznaczają składowe pola prędkości wiatru, czyli w punkcie \( (x,y) \) wiatr wieje w kierunku równoległym do osi \( x \) z prędkością \( \beta_x(x,y) \) a w kierunku równoległym do osi \( y \) z prędkością \( \beta_y(x,y) \).
Z drugiej strony człon \( \epsilon \left(\frac{\partial^2 u }{\partial x^2 }+\frac{\partial^2 u }{\partial y^2 }\right) \) oznacza dyfuzję, czyli roznoszenie zanieczyszczeń na zasadzie dyfuzji.
Jeśli zostanę zamknięty w pokoju z palaczem, to nawet jeśli nie będzieprzeciągu, dym papierosu dotrze do mojego nosa na zasadzie mechanizmu dyfuzji. Za prędkość tej dyfuzji odpowiada współczynnik \( \epsilon \). Jeśli dodatkowo w pokoju było by otwarte okna i drzwi, to dym papierosowy byłby również popychany przez wiatr wiejący z prędkością na przykład \( (1,0) \) (jeśli okno-drzwi leżą w płaszczyźnie równoległe do osi \( x \) i przeciąg wieje z prędkością 1 m/s).
Prawa strona, czyli funkcja \( f(x,y) \) oznacza źródło zanieczyszczeń, czyli jeśli w punkcje \( (x,y) \) znajduje się komin który emituje zanieczyszczenie, wówczas \( f(x,y)>0 \), jeśli natomiast w punkcie \( (x,y) \) znajduje się zasysacz zanieczyszczeń (ostatnio takie zasysasze zanieczyszczeń konstrują Chińczycy w miastach) wówczas \( f(x,y)<0 \). Z koleji \( f(x,y)=0 \) oznacza że w punkcie \( (x,y) \) nie pojawia się ani nie ubywa zanieczyszczeń.
Zauważmy że człon dyfuzji zawiera Laplasjan, a człon adwekcji, iloczyn skalarny wektora \( \beta \) i wektora prędkości wiatru. Możemy więc nasz problem przepisać tak:
\( \beta(x,y) \cdot \nabla u(x,y) - \epsilon \Delta u(x,y) = f(x,y) \)
Sformułowanie słabe uzyskuje się w sposób następujący. Całkujemy i przemnażamy nasze równanie przez wybrane funkcje \( v(x,y) \) zwane funkcjami testującymi których używać będziemy do uśredniania naszego równania w obszarze na którym funkcje te są określone
\( \int_{\Omega} \beta(x,y) \cdot \nabla u(x,y) -\int_{\Omega} \epsilon \Delta u (x,y) v(x,y) dxdy = \int_{\Omega} f(x,y) v(x,y) dxdy \)
W podobny sposób jak miało to miejsce w przypadku problemu projekcji bitmapy, każdy wybór funkcji testującej \( v(x,y) \) służącej do uśredniania naszego problemu w obszarze w którym funkcja testująca jest określona, daje nam jedno równanie. Różne wybory funkcji testującej \( v(x,y) \) dadzą nam więc różne równania .
Nasz problem który chcemy rozwiązać wygląda więc tak: Szukamy koncentracji zanieczyszczeń, funkcji \( \Omega \ni (x,y) \leftarrow u(x,y) \in R \) takiej że
\( \int_{\Omega} \beta_x(x,y) \frac{\partial u(x,y) }{\partial x } dxdy \int_{\Omega } \beta_y(x,y) \frac{\partial u(x,y) }{\partial y } dxdy \\ -\int_{\Omega } \frac{\partial^2 u(x,y) }{\partial x^2 } v(x,y) dxdy -\int_{\Omega } \frac{\partial^2 u(x,y) }{\partial y^2 } v(x,y) dxdy = \\ \int_{\Omega } f(x,y) v(x,y) dxdy \)
dla dowolnych funkcji testujących \( \Omega \ni (x,y) \rightarrow v(x,y) \in R \) (oczywiście nasza koncentracja zanieczyszczeń oraz funkcja testująca muszą być na tyle regularne, żeby dało się policzyć całki).
Oznaczamy nasz problem:
Znajdź \( \Omega \ni (x,y) \rightarrow u(x,y) \in R \) takie że
\( a(u,v)=l(v) \\ a(u,v) = \int_{\Omega} \beta_x(x,y) \frac{\partial u(x,y)}{\partial x} dxdy+ \int_{\Omega } \beta_y(x,y) \frac{\partial u(x,y)}{\partial y } dxdy \\ -\int_{\Omega}\epsilon \frac{\partial^2 u(x,y)}{\partial x^2 } v(x,y) dxdy -\int_{\Omega} \epsilon\frac{\partial^2 u(x,y)}{\partial y^2 } v(x,y) dxdy \\ l(v) = \int_{\Omega} f(x,y) v(x,y) dxdy \)
Względem członu dyfuzji stosujemy wzór na całkowanie przez części, podobnie jak dla równań transportu ciepła. Robimy to dlatego żeby zredukować stopień pochodnych operujących na polu koncentracji zanieczyszczeń.
\( a(u,v) = \int_{\Omega} \beta_x(x,y) \frac{\partial u(x,y)}{\partial x } dxdy + \int_{\Omega} \beta_y(x,y) \frac{\partial u(x,y)}{\partial y } dxdy \\ -\int_{\Omega} \epsilon \frac{\partial^2 u(x,y)}{\partial x^2 } v(x,y) dxdy -\int_{\Omega} \epsilon \frac{\partial^2 u(x,y)}{\partial y^2 } v(x,y) dxdy \\ = \int_{\Omega} \beta_x(x,y) \frac{\partial u(x,y) }{\partial x } dxdy +\int_{\Omega} \beta_y(x,y) \frac{\partial u(x,y)}{\partial y } dxdy \\ +\int_{\Omega} \epsilon \frac{\partial u(x,y) }{\partial x} \frac{\partial v(x,y)}{\partial x } dxdy +\int_{\Omega} \epsilon \frac{\partial u(x,y)}{\partial y } \frac{\partial v(x,y)}{\partial y } dxdy -\int_{\partial \Omega } \frac{\partial u }{\partial n } v dS \)
Zwyczajowo nasz problem adwekcji-dyfuzji, podobnie jak problem transportu ciepła, należy wyposażyć w warunki brzegowe. Dzielimy brzeg \( \partial \Omega = \Gamma_D \cup \Gamma_N \) i wprowadzamy
warunek brzegowy Dirichleta, na \( \Gamma_D \) części brzegu \( \partial \Omega \), który mówi mi o koncentracji zanieczyszczeń na brzegu obszaru \( u(x,y)=p(x,y) \textrm{ dla }(x,y) \in \Gamma_D \)
gdzie \( p(x,y) \) to zadana pochodna koncentracji zanieczyszczeń w kierunku normalnym do brzegu obszaru.
Ponadto wprowadzamy warunek brzegowy Neumanna, który opisuje nam prędkość zmian zanieczyszczeń na \( \Gamma_N \) fragmencie brzegu
\( \frac{\partial u(x,y)}{\partial n}=g(x,y) \textrm{ dla }(x,y) \in \Gamma_N \)
gdzie \( g(x,y) \) to strumień stężenia zanieczyszczeń na brzegu.
Równania konwekcji dyfuzji dla dużej różnicy wartości pomiędzy wektorem adwekcji a współczynnikiem dyfuzji (np. dwa rzędy wielkości lub więcej, czyli np. \( \beta=(1,0), \epsilon=10^{-2} \)) wymagają specjalnych metod stabilizacji opisanych w rozdziale piątym. Dla małych różnic wartości tradycyjna lub izogeometryczna metoda elementów skończonych działa poprawnie.